GetGeoreferenceFromNCdataSet Subroutine

SUBROUTINE GetGeoreferenceFromNCdataSet &
!
(ncId, varId, cellsize, xll, yll, grid_mapping, Iraster, Jraster, point2raster)

USE StringManipulation, ONLY: &
!Imported routines:
StringCompact, &
StringToLower


IMPLICIT NONE

! Arguments with intent (in):
INTEGER (KIND = short), INTENT(IN) :: ncId  !!NetCdf Id for the file
INTEGER (KIND = short), INTENT(IN) :: varId !!variable Id
LOGICAL, INTENT(IN) :: point2raster

! Arguments with intent(out):
REAL (KIND = float), INTENT(OUT) :: cellsize !!cell resolution
REAL (KIND = float), INTENT(OUT) :: xll !!x  coordinate of lower left corner of map
REAL (KIND = float), INTENT(OUT) :: yll !!y  coordinate of lower left corner of map
TYPE (CRS),          INTENT(OUT) :: grid_mapping !!coordinate reference system    

! Arguments with intent(inout): 
INTEGER (KIND = short), INTENT(INOUT) :: Iraster (:,:)
INTEGER (KIND = short), INTENT(INOUT) :: Jraster (:,:)   

! Local variables:
INTEGER (KIND = short) :: x, y !!number of columns and rows
INTEGER (KIND = short) :: idx, idy
REAL (KIND = float), ALLOCATABLE :: xArray (:), yArray (:)
REAL (KIND = float), ALLOCATABLE :: xMatrix (:,:), yMatrix (:,:)
REAL (KIND = float), ALLOCATABLE :: xMatrix2 (:,:), yMatrix2 (:,:)
REAL (KIND = float) :: jdim, idim !! x and y resolution
REAL (KIND = float) :: xmin !! minimum value of x coordinate in xMatrix
REAL (KIND = float) :: xmax !! maximum value of x coordinate in xMatrix
REAL (KIND = float) :: ymin !! minimum value of y coordinate in yMatrix
REAL (KIND = float) :: ymax !! maximum value of y coordinate in yMatrix
INTEGER (KIND = short) :: ncStatus !!error code return by NetCDF routines
INTEGER (KIND = short) :: i, j
INTEGER (KIND = short) :: r, c
CHARACTER (LEN = 80) :: attribute
CHARACTER (LEN = 100)  :: variableName
CHARACTER (LEN = 1) :: shp !!define if x and y are vectors or matrix
INTEGER (KIND = short) :: mappingId !!grid mapping variable Id
INTEGER (KIND = short) :: refSystem !!reference system
INTEGER (KIND = short) :: datum !!datum code
REAL (KIND = float) :: lat0 !!latitude of projection origin
REAL (KIND = float) :: centM !!longitude of central meridian
REAL (KIND = float) :: primeM !!longitude of prime meridian
REAL (KIND = float) :: falseE
REAL (KIND = float) :: falseN
REAL (KIND = float) :: scale_factor
REAL (KIND = float) :: tol = 0.01 !!tolerance for checking regularly spaced grid
INTEGER (KIND = short)  :: latlon
!------------end of declaration------------------------------------------------

!get length of x and y sizes and related id
CALL GetXYSizesFromFile (ncId, x, y, idx = idx, idy = idy, shape = shp, latlon = latlon)


IF (shp == 'v' ) THEN !x and y values are stored in vectors

    !allocate space to read coordinates
    ALLOCATE ( xArray (x) )
    ALLOCATE ( yArray (y) )

    !read x coordinates
    ncStatus = nf90_get_var (ncId, idx, xArray)
    CALL ncErrorHandler (ncStatus)

    !read y coordinates
    ncStatus = nf90_get_var (ncId, idy, yArray)
    CALL ncErrorHandler (ncStatus)
    
    !calculate and check x spatial resolution
    jdim = xArray (2) - xArray (1)
    DO i = 3, x
      IF ( ABS((xArray (i) - xArray (i-1) - jdim ) / jdim) > tol) THEN
        CALL Catch ('error', 'GridLib', 'x dimension not regularly spaced'  )
      END IF
    END DO

    !calculate and check y spatial resolution
    idim = yArray (2) - yArray (1) 
    DO i = 3, y
      IF ( ABS((yArray (i) - yArray (i-1) - idim) / idim) > tol) THEN
        CALL Catch ('error', 'GridLib', 'y dimension not regularly spaced'  )
      END IF
    END DO

    !set cellsize
    IF (ABS((jdim - idim)/jdim) < tol) THEN
      cellsize = (jdim + idim) / 2.
    ELSE
      CALL Catch ('error', 'GridLib', 'x resolution different from y resolution' )
    END IF
     
    !calculate xll and yll
    xll = MINVAL (xArray) - cellsize / 2.
    yll = MINVAL (yArray) - cellsize / 2.

    !deallocate arrays
    DEALLOCATE ( xArray )
    DEALLOCATE ( yArray )
    
ELSE IF (shp == 'm') THEN !x and y values are stored in matrix

    !allocate space to read coordinates
    ALLOCATE ( xMatrix (x,y) )
    ALLOCATE ( yMatrix (x,y) )
    
    ALLOCATE ( xMatrix2 (y,x) )
    ALLOCATE ( yMatrix2 (y,x) )
    
    !read x coordinates
    ncStatus = nf90_get_var (ncId, idx, xMatrix)
    CALL ncErrorHandler (ncStatus)
    
    !read y coordinates
    ncStatus = nf90_get_var (ncId, idy, yMatrix)
    CALL ncErrorHandler (ncStatus)
    
    !calculate x spatial resolution
    xmin = MINVAL(xMatrix)
    xmax = MAXVAL(xMatrix)
    jdim = (xmax - xmin) / (x - 1)
    
    !calculate y spatial resolution
    ymin = MINVAL(yMatrix)
    ymax = MAXVAL(yMatrix)
    idim = (ymax - ymin) / (y - 1)
    
    !set cellsize
    !cellsize = (idim + jdim) / 2.
    cellsize = MAX(idim , jdim)
    
    !calculate xll and yll
    xll = xmin - cellsize / 2.
    yll = ymin - cellsize / 2.
    
    !swap x and y matrix
    CALL SwapGridRealForward (xMatrix, xMatrix2, latlon)
    CALL SwapGridRealForward (yMatrix, yMatrix2, latlon)
    
    !build Iraster and Jraster matrix
    IF (point2raster) THEN
      Iraster = -9999
      Jraster = -9999
      DO i = 1, y
        DO j = 1, x
        
            c = INT ( (xMatrix2(i,j) - xll) / cellsize ) + 1 
            r = INT ( (yll + y * cellsize - yMatrix2(i,j)) / cellsize ) + 1
               
            
             IF (r >= 1 .AND. r <= y .AND. c >= 1 .AND. c <= x ) THEN
             
               Iraster(i,j) = r
               Jraster(i,j) = c
             
             END IF
        
        END DO
      END DO

      
    END IF
    
    !deallocate
    DEALLOCATE (xMatrix)
    DEALLOCATE (yMatrix)
    DEALLOCATE (xMatrix2)
    DEALLOCATE (yMatrix2)
  
END IF



!read coordinate reference system
attribute = ''
ncStatus = nf90_get_att (ncId, varid = varId, name = 'grid_mapping', &
                           values = attribute)

IF (ncStatus == nf90_noerr) THEN !check if grid_mapping exists
  !retrieve varId of grid_mapping variable
  ncStatus = nf90_inq_varid (ncId, attribute, mappingId)
  CALL ncErrorHandler (ncStatus)
  
  !retrieve datum EPSG code (extension to CF conventions)
  ncStatus = nf90_get_att (ncId, varid = mappingId, name = 'datum_code', &
                           values = datum)
  !if datum is not specified, datum_code does not exist, set default to WGS84 
  IF (ncStatus /= nf90_noerr) THEN
    datum = WGS84
  END IF     
  
  !retrieve grid mapping name
  ncStatus = nf90_get_att (ncId, varid = mappingId, name = 'grid_mapping_name', &
                           values = attribute)
  
  IF (ncStatus /= nf90_noerr) THEN !grid_mapping_name not found
      !set default to geodetic WGS84
       CALL SetCRS (GEODETIC, WGS84, grid_mapping)
       CALL SetGeodeticParameters (grid_mapping, prime_meridian = 0.0)
  ELSE
 
      !retrieve reference system parameters                    
      IF ( StringToLower (StringCompact (attribute) ) == 'transverse_mercator' ) THEN
        CALL SetCRS (TM, datum, grid_mapping )
         ncStatus = nf90_get_att (ncId, varid = mappingId, &
                                  name = 'longitude_of_central_meridian', &
                                  values = centM)
         ncStatus = nf90_get_att (ncId, varid = mappingId, &
                                  name = 'central_meridian', &
                                  values = centM)  
         ncStatus = nf90_get_att (ncId, varid = mappingId, &
                                  name = 'latitude_of_projection_origin', &
                                  values = lat0)
         ncStatus = nf90_get_att (ncId, varid = mappingId, &
                                  name = 'scale_factor_at_central_meridian', &
                                  values = scale_factor)
          ncStatus = nf90_get_att (ncId, varid = mappingId, &
                                  name = 'scale_factor', &
                                  values = scale_factor)
         ncStatus = nf90_get_att (ncId, varid = mappingId, &
                                  name = 'false_easting', &
                                  values = falseE)
         ncStatus = nf90_get_att (ncId, varid = mappingId, &
                                  name = 'false_northing', &
                                  values = falseN)
         CALL SetTransverseMercatorParameters (grid_mapping, lat0, centM, &
                                               falseE, falseN, scale_factor)
      ELSE IF ( StringCompact (attribute) == 'latitude_longitude' ) THEN
         CALL SetCRS (GEODETIC, datum, grid_mapping )
         ncStatus = nf90_get_att (ncId, varid = mappingId, &
                                  name = 'longitude_of_prime_meridian', &
                                  values = primeM)
     
         CALL SetGeodeticParameters (grid_mapping, primeM)
      ELSE
      !other reference systems not yet supported
      END IF
   END IF 
ELSE !CRS not specified: set default to geodetic WGS84
  !retrieve name of variable
   ncStatus = nf90_inquire_variable(ncId, varId, name = variableName)
   CALL ncErrorHandler (ncStatus)
   CALL Catch ('warning', 'GridLib',        &
    'grid_mapping not found in variable: ',  &
    argument = variableName ) 
    
    CALL Catch ('info', 'GridLib',        &
    'set default coordinate reference system to geodetic wgs84') 
    
    CALL SetCRS (GEODETIC, WGS84, grid_mapping)
    CALL SetGeodeticParameters (grid_mapping, prime_meridian = 0.0)

END IF                           

RETURN
END SUBROUTINE GetGeoreferenceFromNCdataSet